MODULE SRFSN_RSN_MOD
CONTAINS
SUBROUTINE SRFSN_RSN(KIDIA,KFDIA,KLON,KLEVSN,PTMST,LLNOSNOW,PCIL,LDLAND,&
                    &PFRSN,PRSNM1M,PSSNM1M,PTSNM1M,PWSNM1M,PWSN,&
                    &PSNOWF,PUSRF,PVSRF,PTSRF,&
                    &YDSOIL,YDCST,PRSN,PDHTSS)


USE PARKIND1 , ONLY : JPIM, JPRB
USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK

USE YOS_SOIL , ONLY : TSOIL 
USE YOS_CST  , ONLY : TCST

! (C) Copyright 2015- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!**** *SRFSN_RSN* - Snow density
!     PURPOSE.
!     --------
!          THIS ROUTINE COMPUTES THE SNOW DENSITY
!          As a single prognostic snowpack for seasonal snow and land ice is used,
!          the two different contributions are weighted by PCIL for sub-grid ice.

!**   INTERFACE.
!     ----------
!          *SRFSN_RSN* IS CALLED FROM *SRFSN_DRIVER*.

!     PARAMETER   DESCRIPTION                                    UNITS
!     ---------   -----------                                    -----

!     INPUT PARAMETERS (INTEGER):
!    *KIDIA*      START POINT
!    *KFDIA*      END POINT
!    *KLON*       NUMBER OF GRID POINTS PER PACKET
!    *KLEVSN*     VERTICAL SNOW LAYERS

!     INPUT PARAMETERS (REAL):
!    *PTMST*      TIME STEP                                      S

!     INPUT PARAMETERS (LOGICAL):
!    *LLNOSNOW*   NO-SNOW/SNOW MASK (TRUE IF NO-SNOW)

!     INPUT PARAMETERS AT T-1 OR CONSTANT IN TIME (REAL):
!    *PSSNM1M*    TOTAL SNOW MASS IN EACH LAYER (per unit area) kg/m**2
!    *PWSNM1M*    LIQUID WATER CONTENT IN SNOW                 kg/m**2
!    *PRSNM1M*    SNOW DENSITY in each layer                   kg/m**3
!    *PTSNM1M*    TEMPERATURE OF SNOW LAYER                    K
!    *PSNOWF*     TOTAL SNOW FLUX AT THE SURFACE              KG/M**2/S
!    *PUSRF*      U-WIND COMPONENT LOWER MODEL LEVEL           m/s
!    *PVSRF*     V-WIND COMPONENT LOWER MODEL LEVEL            m/s
!    *PTSRF*      TEMPERATURE LOWER MODEL LEVEL                 K

!     INPUT PARAMETERS AT T
!    *PWSN* LIQUID WATER CONTENT IN SNOW                 kg/m**2

!     OUTPUT PARAMETERS AT T+1 (UNFILTERED,REAL):
!    *PRSN*       SNOW DENSITY in each layer                   kg/m**3

!     METHOD.
!     -------
!          

!     EXTERNALS.
!     ----------
!          NONE.

!     REFERENCE.
!     ----------
!          

!     Modifications:
!     Original   E. Dutra      ECMWF     04/12/2015
!                G. Arduini    ECMWF     01/09/2021
!     Modified   G. Arduini    ECMWF     Sept 2024    snow over land-ice
!     ------------------------------------------------------------------

IMPLICIT NONE

! Declaration of arguments 
INTEGER(KIND=JPIM), INTENT(IN)   :: KIDIA
INTEGER(KIND=JPIM), INTENT(IN)   :: KFDIA
INTEGER(KIND=JPIM), INTENT(IN)   :: KLON
INTEGER(KIND=JPIM), INTENT(IN)   :: KLEVSN
REAL(KIND=JPRB)   , INTENT(IN)   :: PTMST
REAL(KIND=JPRB)   , INTENT(IN)   :: PCIL(:)
LOGICAL           , INTENT(IN)   :: LLNOSNOW(:) 
LOGICAL           , INTENT(IN)   :: LDLAND(:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PFRSN(:)

REAL(KIND=JPRB)   , INTENT(IN)   :: PRSNM1M(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PSSNM1M(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PTSNM1M(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PWSNM1M(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PWSN(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PSNOWF(:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PUSRF(:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PVSRF(:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PTSRF(:)
TYPE(TSOIL)       , INTENT(IN)   :: YDSOIL
TYPE(TCST)        , INTENT(IN)   :: YDCST

REAL(KIND=JPRB)   , INTENT(OUT)  :: PRSN(:,:)

REAL(KIND=JPRB)   , INTENT(OUT)  :: PDHTSS(:,:,:)

! Local variables 
INTEGER(KIND=JPIM)  :: JL,JK

REAL(KIND=JPRB)     :: ZRHOFRESH   ! FRESH SNOW DENSITY
REAL(KIND=JPRB)     :: ZRSTAR(KLEVSN)      ! NEW DENSITY AFTER ACCOUNTING FOR SNOWFALL 
REAL(KIND=JPRB)     :: ZSNOWFM     ! SNOWFALL MASS 
REAL(KIND=JPRB)     :: ZRSNDTOVER  !  OVERBURDEN
REAL(KIND=JPRB)     :: ZRSNDTDEST  ! METAMORPHISM
REAL(KIND=JPRB)     :: ZRSNDTWINT  ! LIQUID WATER INTERCEPTION
REAL(KIND=JPRB)     :: ZRSNDTDESTCV
REAL(KIND=JPRB)     :: ZEPS,ZADD,ZTMST

REAL(KIND=JPRB)     :: ZRSNDTWINDCOMP ! WIND-DRIVEN COMPACTION
REAL(KIND=JPRB)     :: ZFWCOMP1, ZFWCOMP2, ZFTAUW, ZTAUW
REAL(KIND=JPRB)     :: ZGAMMAMOB(KLEVSN), ZGAMMAW(KLEVSN), ZDSN(KLEVSN)
REAL(KIND=JPRB)     :: ZSD_THR, ZDSN_SUM
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SRFSN_RSNL_MOD:SRFSN_RSN',0,ZHOOK_HANDLE)

!    -----------------------------------------------------------------
ASSOCIATE(RHOMINSN=>YDSOIL%RHOMINSN, RHOMINSNA=>YDSOIL%RHOMINSNA, &
 & RHOMINSNB=>YDSOIL%RHOMINSNB, RHOMINSNC=>YDSOIL%RHOMINSNC, &
 & RHOMINSND=>YDSOIL%RHOMINSND, RSNDTDESTA=>YDSOIL%RSNDTDESTA,&
 & RSNDTDESTB=>YDSOIL%RSNDTDESTB, RSNDTDESTC=>YDSOIL%RSNDTDESTC, &
 & RSNDTDESTC_ML=>YDSOIL%RSNDTDESTC_ML, &
 & RSNDTDESTROI=>YDSOIL%RSNDTDESTROI, RSNDTOVERA=>YDSOIL%RSNDTOVERA, &
 & RSNDTOVERB=>YDSOIL%RSNDTOVERB, RSNDTOVERC=>YDSOIL%RSNDTOVERC ,&
 & RG=>YDCST%RG,RTT=>YDCST%RTT, RSNDAMOB=>YDSOIL%RSNDAMOB,       &
 & RSNDMOB=>YDSOIL%RSNDMOB, RSNDAW=>YDSOIL%RSNDAW,               &
 & RSNDBW=>YDSOIL%RSNDBW, RSNDKV=>YDSOIL%RSNDKV,                 &
 & RSNDATAUW=>YDSOIL%RSNDATAUW, RSNDBTAUW=>YDSOIL%RSNDBTAUW,     &
 & RSNDWCOMPMAX=>YDSOIL%RSNDWCOMPMAX,RHOMAXSN_NEW=>YDSOIL%RHOMAXSN_NEW )

ZTMST = 1.0_JPRB/PTMST 

ZEPS=10._JPRB*10._JPRB**(-PRECISION(1._JPRB))
DO JL=KIDIA,KFDIA
  IF (LLNOSNOW(JL)) THEN
    PRSN(JL,:)=RHOMINSN
  ELSE
    ZRSTAR(:) = PRSNM1M(JL,:)
    
  !*  FRESH SNOW DENSITY
    ZRHOFRESH = RHOMINSNA + RHOMINSNB*(PTSRF(JL)-RTT) + &
        & RHOMINSNC*(PUSRF(JL)**2+PVSRF(JL)**2)**0.25_JPRB
    ZRHOFRESH=MAX(ZRHOFRESH,RHOMINSND)
    ! Assuming fresh snow over full glacier points at 300 kg m-3
    IF (PCIL(JL)>1._JPRB-ZEPS .AND. LDLAND(JL))THEN
      !*ZRHOFRESH=300._JPRB !Alexander et al. JGR 2019
      ZRHOFRESH=280._JPRB !Alexander et al. JGR 2019
    ENDIF
  !* FIRST GUESS SNOW DENSITY AFTER SNOWFALL
    ZSNOWFM = PTMST*PSNOWF(JL) 
    ZRSTAR(1) = ( PSSNM1M(JL,1) + ZSNOWFM ) / MAX( ZEPS,&
    & ( ZSNOWFM/ZRHOFRESH + PSSNM1M(JL,1)/PRSNM1M(JL,1) ) )
    
    ZGAMMAMOB(1:KLEVSN)=0._JPRB
    ZGAMMAW(1:KLEVSN)=0._JPRB
    ZDSN(1:KLEVSN)=0._JPRB
  !* Compute mobility index and wind-blowing index for computation of tau in wind compaction
    DO JK=1,KLEVSN
      ZGAMMAMOB(JK)=RSNDAMOB*( 1.0_JPRB - MAX( 0._JPRB, (PRSNM1M(JL,JK) - RHOMINSND)/RSNDMOB ) )
      ZGAMMAW(JK)=1._JPRB - RSNDAW*MAX(ZEPS,EXP(-RSNDBW*SQRT(PUSRF(JL)**2+PVSRF(JL)**2))) + ZGAMMAMOB(JK)
      ZDSN(JK)=PSSNM1M(JL,JK) / MAX(ZEPS,PRSNM1M(JL,JK)*PFRSN(JL))
    ENDDO
    ! COMPACTION RATES PER LAYER
    ZADD=0._JPRB      
    ZFWCOMP1=0._JPRB
    ZFWCOMP2=0._JPRB
    ZTAUW=0._JPRB
    DO JK=1,KLEVSN
      IF (PSSNM1M(JL,JK) > ZEPS ) THEN
       !*    OVERBURDEN
        ZRSNDTOVER=(ZADD+0.5*PSSNM1M(JL,JK))*RG/(RSNDTOVERA*&
        & EXP(RSNDTOVERB*(RTT-PTSNM1M(JL,JK))+RSNDTOVERC*PRSNM1M(JL,JK)))

       !*    METAMORPHISM  
        ZRSNDTDEST=RSNDTDESTA*EXP(-RSNDTDESTB*(RTT-PTSNM1M(JL,JK))-&
      & RSNDTDESTC_ML(JK)*MAX(0.0_JPRB,PRSNM1M(JL,JK)-RSNDTDESTROI))

        !*    COMPACTION DUE TO LIQUID WATER INTERCEPTION (OR MELTING RETAINED)
        ZRSNDTWINT = ZTMST*MAX(0._JPRB,PWSN(JL,JK)-PWSNM1M(JL,JK))/MAX(ZEPS,PSSNM1M(JL,JK))

        !* COMPACTION DUE TO WIND-INDUCED DENSIFICATION
        ZFWCOMP1=MAX( 0._JPRB, ZGAMMAW(JK))
        ZFWCOMP2=ZFWCOMP2 + (ZDSN(JK)*( RSNDBTAUW - ZGAMMAW(JK) ))
        ZFTAUW=ZFWCOMP1 * MAX(ZEPS,EXP(-RSNDATAUW * ZFWCOMP2))
        ZTAUW=2._JPRB*86400.0_JPRB/MAX(ZEPS,ZFTAUW)
        
        ZRSNDTWINDCOMP=MAX(0._JPRB, (RSNDWCOMPMAX-PRSNM1M(JL,JK))/ZTAUW)

        !* Final computation            
        ZRSTAR(JK) = ZRSTAR(JK)+ZRSTAR(JK)*(ZRSNDTOVER+ZRSNDTDEST+ZRSNDTWINT)*PTMST + ZRSNDTWINDCOMP*PTMST
        ZRSTAR(JK) = MIN(500._JPRB,ZRSTAR(JK))
        ZRSTAR(JK) = MAX(RHOMINSND,ZRSTAR(JK))

        ZADD = ZADD + PSSNM1M(JL,JK)
      ENDIF
    ENDDO
    IF (PCIL(JL)<1._JPRB-ZEPS .AND. LDLAND(JL))THEN
      PRSN(JL,:)=MIN(RHOMAXSN_NEW,ZRSTAR(:))
      PRSN(JL,:)=MAX(RHOMINSND,ZRSTAR(:))
    ELSE
      PRSN(JL,1)=MIN(315._JPRB,ZRSTAR(1))  !Fausto et al. Frontiers 2018 ~(315-44)
      ZSD_THR=1._JPRB
      ZDSN_SUM=PSSNM1M(JL,1) / PRSNM1M(JL,1)
      DO JK=2,KLEVSN
            PRSN(JL,JK)=MAX(PRSN(JL,1), MIN(RHOMAXSN_NEW, ZRSTAR(JK)) ) ! layers below 1m are evolving
      ENDDO
    ENDIF


  ENDIF
ENDDO
  
END ASSOCIATE
!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SRFSN_RSNL_MOD:SRFSN_RSN',1,ZHOOK_HANDLE)

END SUBROUTINE SRFSN_RSN
END MODULE SRFSN_RSN_MOD
